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Abstract 

In this paper, we solve the field equations in metric f(R) gravity 
for Bianchi type VIo spacetime and discuss evolution of the expand- 
ing universe. We find two types of non-vacuum solutions by taking 
isotropic and anisotropic fluids as the source of matter and dark en- 
ergy. The physical behavior of these solutions is analyzed and com- 
pared in the future evolution with the help of some physical and geo- 
metrical parameters. It is concluded that in the presence of isotropic 
fluid, the model has singularity at t = and represents continuously 
expanding shearing universe currently entering into phantom phase. 
In anisotropic fluid, the model has no initial singularity and exhibits 
the uniform accelerating expansion. However, the spacetime does not 
achieve isotropy as t — > oo in both of these solutions. 
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1 Introduction 

Extended theories of gravity have become a paradigm in bringing suitable 
cosmological models where a late time accelerated expansion can be achieved. 
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The paradigm, consisting of higher order curvature, establishes the f(R) the- 
ory of gravity. This theory has many applications in cosmology and gravity 
such as inflation, dark energy (DE), local gravity constraints, cosmological 
perturbations and spherically symmetric solutions in weak and strong grav- 
itational backgrounds. The main motivation of this theory comes from the 
fact that every unification of fundamental interaction exhibits effective ac- 
tions containing higher order terms in the curvature invariants. This strategy 
was adopted in the study of quantum field theory in curved spacetimes [1] 
as well as in the Lagrangian of string and Kaluza-Klein theories [2] . 

Higher order terms always give an even number as an order of the field 
equations. For example, R 2 term produces fourth order field equations [3J, 
term ROR (where □ = V^V^) gives sixth order field equations jU simi- 
larly, R\3 2 R yields eights order field equations [6] and so on. Using conformal 
transformation, the term with second derivative corresponds to a scalar field. 
For instance, fourth order gravitational theory corresponds to Einstein the- 
ory with one scalar field, sixth order gravity corresponds to Einstein gravity 
with two scalar fields, etc. [HE]. In this context, it is easy to show that f(R) 
gravity is equivalent to scalar tensor theory as well as to general relativity 
(GR) with an ideal fluid [8]. 

In an isotropic and homogeneous spacetimes, the Einstein field equations 
give rise to the Friedmann equations which are used to describe the evo- 
lution of the universe. However, the speedy development in observational 
cosmology shows that the universe has undergone two phases of cosmic ac- 
celeration. The first phase is called inflation [9]- [12] which is believed to have 
occurred earlier than matter domination [13]- [15]. This accelerating phase is 
expected to solve horizon problems involved in the big bang cosmology and 
to explain nearly flat spectrum of temperature anisotropies observed in Cos- 
mic Microwave Background (CMB) radiations [TB]. The second accelerating 
phase has occurred after the matter domination. 

The first model of inflation with R + aR 2 , proposed by Starobinsky [9], 
can lead to an accelerated expansion of the universe due to aR 2 term. The 
unknown component giving rise to this late-time acceleration is called dark 
energy (DE) [T7]-[2T]. The cosmic acceleration has been confirmed by many 
observational sources such as Supernovae la (SNela) [22]- [21], Large scale 
structure (LSS) [23 [26] , CMB (2Z]-[29J and baryon acoustic oscillations 
(BAO) [301 13T] . The discovery of DE stimulated the idea that cosmic ac- 
celeration might originate from some extension of GR. Dark energy models 
based on f(R) theory have been extensively studied to explain the late-time 
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acceleration. 

Most of the exact solutions in this theory have been discussed in spherical 
symmetry. Multamdki and Vilja [32j studied static spherically symmetric 
vacuum as well as non- vacuum solutions by taking perfect fluid (33] . Carames 
and Bezerra [31] found spherically symmetric vacuum solutions in higher 
dimensions. Capozziello et al. [35J analyzed spherically symmetric solutions 
using Noether symmetry. The structure of relativistic stars in f(R) theory 
has been discussed by many authors [36]- [30]. In a recent paper, we have 
found static spherically symmetric solution in the presence of dust matter and 
discussed energy-momentum distribution for some well-known f(R) models 
[4"T] using Landau-Lifshitz complex. We have also checked the stability and 
constant curvature conditions of these models. 

In cylindrical symmetry, Azadi et al. [42] have studied solutions in Weyl 
coordinates. They have shown that constant curvature solutions reduce to 
only one member of the Tian family in GR. Momeni (33] has found that exact 
constant scalar curvature solution in cylindrical symmetry is applicable to 
the exterior of a string. Sharif and his collaborators [331 132] have investigated 
the exact solutions of plane Bianchi models for both vacuum and non-vacuum 
cases. Hollestein and Lobo [46J explored exact solutions of the field equations 
coupled to non-linear electrodynamics. 

It has been observed that some large-angle anomalies appear in CMB ra- 
diations which violate the statistical isotropy of the universe [4T] 145] . Plane 
Bianchi models (which are homogeneous but not necessarily isotropic) seem 
to be the most promising explanation of these anomalies. Jaffe et al. [49]- [51] 
investigated that removing a Bianchi component from the WMAP data can 
account for several large-angle anomalies leaving the universe to be isotropic. 
Thus the universe may have achieved a slight anisotropic geometry in cos- 
mological models regardless of the inflation. Further, these models can be 
classified according to whether anisotopy occurs at an early stage or at later 
times of the universe. The models for the early stage can be modified in a 
way to end inflation with a slight anisotropic geometry [52]. For the later 
class, the isotropy of the universe, achieved during inflation, can be distorted 
by modifying DE [23] 123] • 

The objective of this paper is to find non-vacuum exact solutions of 
Bianchi type VIq model in the metric f{R) gravity. The organization of 
the paper is as follows. In section 2, we present the field equations and some 
dynamical quantities describing the evolution of the universe. Section 3 pro- 
vides solution of the field equations in the presence of perfect fluid. In section 
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4, we find solution for anisotropic fluid. In the last section 5, summary and 
comparison of both the solutions is given. 



2 The Model and the Field Equations 

The line element for homogeneous and anisotropic Bianchi type VIq space- 
time is given by 

ds 2 = dt 2 - A 2 (t)dx 2 - e 2ax B 2 (t)dy 2 - e~ 2ax C 2 {t)dz 2 , (2.1) 

where the scale factors A, B and C are functions of cosmic time t only and a 
is a non-zero constant. There are two formalisms which are applied to derive 
the field equations in f(R) gravity. One is the standard metric formalism 
while another is Palatini formalism. Most of the work in this theory has been 
done by using the former formalism, where the action with S^ m ' as a matter 
part, is given as follows [55] 

S = YkS d4x ^f{R) + S^ m) . (2.1) 

Its variation with respect to the metric tensor yields the following set of field 
equations 

F{R)Rp, - \f(R)g^ - V M V„F(i2) + g^DF(R) = kT^, (2.2) 

where F(R) = df(R)/dR and f(R) is a function of the Ricci scalar and 
describes all kinds of matter including non-relativistic (cold) dark matter. 
Taking trace of the above equation (with k — 1), we obtain f(R) 

f[R) = mR+^m-r (2 3) 

The scalar curvature for Bianchi type VIq model is given by 



R 



A B C a 2 AB AC BC 
-r + TT + ^- -TTT + -r^ + 



ABC A 2 AB AC BC 



(2.4) 



The directional Hubble parameters along x, y and z directions, the aver- 
age scale factor a, the deceleration parameter q, the expansion scalar B, the 
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shear scalar a and the mean Hubble parameter H are respectively given as 

(ABC) l/ \ (2.5) 
o 2 = \o ab a a \ (2.6) 
(2.7) 

Here u a is the four velocity vector and a ab is the shear tensor. Berman [561157] 
introduced the variation of mean Hubble parameter as follows 

H = ka~ n , (2.8) 

where k > and n ^ 0. Inserting this value in Eq.f l2.6p . we obtain 

a = de kt , q=~l for n = 0, (2.9) 

a = (nkt + c 2 ) 1/n , q = n-l for n ^ 0, (2.10) 

where C\ and c 2 are positive constants of integration. For n < 1, these 
equations represent an accelerated expansion of the universe with a — > oo 
as t — )• oo and supports the observations of the Type la supernova [22} [23] 
and WMAP data [271 128] . To examine whether expansion of the universe is 
anisotropic or not, we define anisotropic expansion parameter as 

If A = 0, then the expansion of the universe is isotropic. Further, any 
anisotropic model of the universe with diagonal energy-momentum tensor 
approaches to isotropy if A — >■ 0, V — > +oo and p > as t — > +oo [SHI 159] - 

3 Solution with Isotropic Fluid 

In this section, we take isotropic perfect fluid whose energy-momentum tensor 
is given by 

T u = diag[p, -p, -p, -p], (3.1) 
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where p is the density and p is the pressure. The corresponding field equations 
become 

A 2a 2 AB AC s „ 1, /t , n « ,S C\- , , 

)F+ 1 -f(R)-F-(j + ^)F=p, (3.4) 

)^+^(^)-^-(^ + |)^=P, (3-5) 

«(§-§) F = °- ( 3 - 6 ) 

The solution of Eq. 03.61) yields 

C = c 3 B, (3.7) 

where C3 > is another constant of integration. Without any loss of gener- 
ality, we take C3 = 1 for the sake of simplicity. Using this value of C in the 
above equations, we obtain 

{ * + § )F+ l m) - { A + ™ } P = (3.8) 
,A 2a 2 2AB,„ 1 ~ 25- , , 

<i-^ + ^r) F+ 2 /(fl) - F -ir F = p - (M) 

Subtracting Eq. fl3.9j) from (13. 8p . Eq. fl3.10j ) from (13.81) respectively and divid- 
ing the resulting equations by F, we have 

-B 2a 2 2AB F AF p + p 

~B~~^ + ^B~F + AF ~ "p"' (3 ' U) 
A B AB B 2 F BF p + p 

~A~B + AB + B 2 ~F + BF ~ ~F~' 

We assume that the scalar expansion is proportional to the shear scalar 
[60]- (62], i.e., O oc a, which leads to a relation between the metric functions 
as follows 

B = A n , (3.13) 
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where n ^ 1 is a positive constant. Using this relation in Eqs.f l3.lTj) and 
fl3TT2|) . it follows that 



A A 2a 2 F AF 

2n(2 - n)— - 2n 1 

V 'A 2 A A 2 F AF 

A 2 . 'A F AF 

2n Al- {n + 1) A-F +n AF 

Subtraction of Eq.f l3.15j) from (13. 14j) yields 



P + P 
F ' 

P + P 



(3.14) 
(3.15) 



A A 2 AF 2a 2 

l_ 2tj I 

A A 2 AF {l-n)A 2 



0. 



We solve this equation using power law relation between F and a 

F = la m , 

where I is the constant of proportionality and m is any real number, 
tuting the value of a, F turns out to be 

F = lA^ m . 

Inserting this value in Eq. (13. 16j) . it follows that 

A 6n + (2n + l)m A 2 2a 2 



A 3 A 2 {l-n)A 2 

Integrating this equation, we obtain 



0. 



(3.16) 

j 

(3.17) 
Substi- 

(3.18) 
(3.19) 



A 



6a 2 



C4 



(1 - Tl) (m + 6n + 2mn) ^2(m+6n+2mn)/3 ' 



(3.20) 



where C4 is another integration constant. With the help of Eqs. fl3.7p and 
(13. 13 j) . the line element (12 .11) reduces to 



ds 1 



(1 - n)(m + 6n + 2mn)A 2 ^ m+&n+2mn ^ z 



g Q ,2 j 42(m+6n+2 m n)/3 + c ^ m + Q n + 2 mn ) 

- e 2ax A 2n dy 2 - c\A 2n e- 2ax dz 2 



dA A - A 2 dx 2 



(3.21) 
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Taking the coordinate transformation (A = t), Eq.( l3.21j) becomes 



(1 — n iim -+- on + zmn ir~ v '" 1 — ~ 9 ~ 9 , 9 

i-i i djf- — t ax 

Q a 2 t 2(m+6n+2mn)/3 + c ^ m + Q n + 2mn) J 

2ax P n dy 2 - c 2 3 P n e- 2ax dz 2 . (3.22) 



n){m + 6n + 2mn)f 2(m+6n+2 ™ )/3 



rfs 2 = 
— e 

Using Eq. @rrgp in Eqs. flBTTIj) and (j3TT5|l . we obtain 

2(2n + l)m (2n + l) 2 m 2 i 2 



[2n(2 - n) + 



3 9 M 2 

(2n + l)mv4 2a 2 p + p 

3 J I~"Z? 



-[2™+- ~ — = . . (3-23) 



-K» + l) + ^=]4 = 7^- (324) 



Since the right hand sides of these equations are same with two unknown 
functions p and p, thus we get a dependent solution. Using equation of state 
(EoS) p = up, where u is the EoS parameter and adding Eqs. (13.231) and 
(I3.24p . the energy density is found to be 

~ (2n + l)m „ „ 

It a z —6a c 4 , 

9 ~ 2(1 + U) ^ (n - l)(m + 6n + 2mn) + £2(m+6n+2mn)/3 

x {6n 2 + m(-2n 2 + — + -) + ~(2n + l) 2 m 2 } 

12 3 9 

9 ,71 + 1 ,271 + 1,771,, 

+ 4Q '(— + <— >3>1- (3 ' 25) 

The positive values of u characterize different kinds of fluids while in order 
to explain DE of the current universe with the help of Bianchi models, one 
can use u with negative values. When w — — 1, the universe passes through 
ACDM epoch. If w < —1, then we live in the phantom-dominated universe 
and for w > —1, the quintessence dark era occurs. Here we analyze the 
behavior of energy density in phantom and quintessence regions for different 
positive values of m as shown in Figure 1 while u — — 1 yields singular 
expression. It can be observed that our solution favors only the phantom 
region because p is negative in the other region which is not feasible. Thus 
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Phantom region 



Quintessence region 



P Behavior of energy density 



P Behavior of energy density 





Figure 1: Behavior of energy density versus cosmic time. The three colored 
lines are sketched for m = 2, m = 3 and m = 4 respectively. 



the energy density will increase and the universe enters into the phantom 
phase. It is worthwhile to mention here that observational analysis of a recent 
supernova strongly supports w < — 1 being EoS parameter for phantom DE 
[62]-[65]. 

For this model, the physical quantities will become 
(2n + 1) 



H = 3t -m+6n+2 m n+i V -2a 2 t 2 (">+ 6 "+ 2m ") + c 4 (n - l)(m + 6n + 2mn), 
V = P n+1 . (3.26) 

The expansion scalar and the shear scalar a become 

(2n + i)(_6a 2 f 2 ( m+6n+2mn )/ 3 + c 4 {n - l)(m + 6n + 2mn)) 1 / 2 

" ~ rn + 6n + 2mn i -i 5 \" - ^'/ 

r 3 ^ 1 

(1 - n)(-6a 2 t 2 ( ,n+6n+2mn ) + c 4 (n - l)(m + 6n + 2mn)) 1 / 2 

°" ~ m+6n + 2mn , \ • (3.28) 

V3 t 3 +1 

We see that the scale factors and volume of the universe are zero at initial 
epoch which shows that the model has point type singularity while they 
continue to increase with time. The Hubble parameter and the expansion 
scalar indicate that expansion rate was rapid at initial times of the big bang 
but it slows down with the passage of time and tends to zero as t — > oo. 
The ratio g indicates that the universe does not achieve isotropy and hence 
model represents continuously expanding, shearing universe from the start 
of the big bang. 
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The scalar curvature and f(R) function turn out to be 

2a 2 (5n +1) , , N9 , 

R = ^-7 - 1 + 2(5n 2 + 6n + (2n + 1) 2 to) 

£ 2 (n-l) 

-2a 2 f 2 ( m+6n+2mn ) + c 4 (n - 1)(to + 6n + 2mn) 
(n - l)(m + 6n + 2mn)f 2 (™+ 6n + 2m "+ 1 ) ' 



= lP n+1)/3 R + 3(n + i; 



—2a 2 c 4 



_ J ^ lV2 (o oq\ 



(n - 1) (to + 6n + 2mn)t 2 t 2 ( m+6n+2mn+1 ^ 

We know that if /(i?) is replaced by i? + A, then f(R) theory corresponds 
to GR, where A is interpreted as the energy density of the vacuum [66J- 
[69] causing expansion in the universe. Thus f(R) may be used to explain 
the present cosmological expansion. We note that the above function shows 
initial time singularity and continuously expanding to infinity. 



4 Solution with Anisotropic Fluid 

Here, we obtain solution of the field equations for anisotropic fluid 

T£ = diag[p, -p x , -p y , -p z ] = diag[l, -u x , -u y , -u z }p, (4.1) 

where p x , p y and p z are pressures and u x , uj y and u z are directional EoS 
parameters on x, y and z axes respectively. We take 

u x = u + 5, Uy = u and co z = u + 7, 

where u is the deviation-free EoS parameter and 5 and 7 (called skewness 
parameters) are deviations from u> on x and z axes respectively. The energy- 
momentum tensor takes the form 

T; = diag[l, -(u + 5), -u>, -(u + 7 )]p. (4.2) 
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Notice that u, 5 and 7 are functions of cosmic time t. The field equations 
will become 

«(§-|)F = 0. (4.7) 

Using solution of Eq. (I4.7j) in the above system and then subtracting 
Eq. fl4.5p from (14. 6ft . we obtain 7 = 0. This indicates that the directional 
EoS parameters u y , u z along y and z axes become equal and hence also 
pressure. Consequently, the field equations turn out to be 

,A 2B, 1 „, N .A 2B ■ . , 

( I + 5" )F+ 2 /(jR)_( I + = (4 ' 8) 

Subtracting Eq. fl4.9|) from (I4.10p and integrating the resulting equation, we 
have 

n* - H v = Yf + Vf I { ^ + (411) 

where C5 is another positive constant of integration. Using Eqs.( l2.1iT) and 
(ICTj) . it follows that 

2 r a 2 F 

A = ^[c 5 + y {-^ + 5p)Vdt?V-*F-\ (4.12) 

If we take 5 = and F(R) = 1, the anisotropy parameter of expansion 
reduces to GR for an isotropic fluid. To avoid the integral term in Eq. 04.12]) . 
we can choose 5 such that 

2a 2 F 

6 = -^. (4.13) 
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The energy- momentum tensor, the anisotropy parameter and Eq. fl4.lip will 
take the form 



2a 2 F 
pA 2 



T; = diag[l, -u + — -u, -u]p, (4.14) 



2 cl 



A = (4 - 15) 



H x ~H y = ^. (4.16) 



Using Eq.f l4.13p in (14. 9p . we obtain the following set of field equation 

(l + f)^>-4 + f>^ = "A (4-17) 
+ ^ + = (4.18) 

( B + iB + iF» F+ 5 /(fl) - F - ( T + ir) F = ^ < 4 - 19 > 

We can solve this system of equations by using exponential and power 
law volumetric expansions corresponding to Eqs.f l2.9p and f l2.10p respectively. 
The exponential expansion (also known as de Sitter expansion) has been 
discussed in detail [70J. Using Eq. fl2.9p . the volume of the universe and F 
can be written as follows 

V = c?e 3fct , (4.20) 
F = lc™e mkt . (4.21) 

Inserting Eqs. fT4~20l and ( OB in f l4~T6l) . we have 

r p — (3+m)kt 

H x ~H y = 5 ic , +m . (4.22) 

Solving the system of equations fl4.17p - fl4.19l) along with Eq.f l4.22p . we obtain 
the scale factors as follows 

kt 2 c 5 e -(3+m)fct 

4 = c 6 e 3 «(3+ m) * c ?+ m } ^ 4 _23) 

B = C = c 7 e ^+™)*4 +m ^ ^ 4 _24) 

where C6 and C7 are positive constants of integration. When m > —3, 
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Behavior for m<-3 




0T5 1 TTs 2 2T5 3 

Figure 2: Behavior of scale factors versus time 



the scale factors admit constant values at initial time afterwards they start 
increasing with cosmic time without any type of initial singularity and finally 
diverges to 00 as t — > 00 (Figure 2a). This shows that at the initial epoch, 
the universe starts with some constant volume and expands exponentially 
approaching to infinite volume. However, for m < —3, the scale factor A 
increases rapidly with time (Figure 2b) while B first increases for a finite 
value of time and then decreases approaching to zero at later times (Figure 
2c). Moreover, the expansion scalar is constant showing that the universe is 
expanding uniformly from t — to t — 00. 

The directional and mean Hubble parameters become 

H x = k + ^-e~^ k \ H v = H z = k %-e~^ kt , 

3l4 +m y 3lcl +m 

H = k. (4.25) 
The anisotropy parameters turns out to be 

A = §L^ 2 ( m+3 ) fci . (4.26) 
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We see that the mean Hubble parameter is constant whereas others are time 
dependent. As time approaches from zero to infinity (for m > —3), the 
directional Hubble parameters are reduced to the mean value. For m < —3, 
parameter along x-axis will increase from the mean Hubble parameter by 
twice a constant factor (in the coefficient of exponential function) whereas 
parameters along y and z-axes decrease by the same factor. Notice that at 
t — 0, the anisotropy parameter measures a constant value while it vanishes 
for m > —3 at infinite time of the universe. This indicates that the universe 
expands isotropically at later times without taking any effect from anisotropy 
of the fluid. However, for m < —3, the anisotropy in the expansion will 
increase with time. 

Using values of the scale factors in Eq. fl4.8p . we obtain the following 
energy density of anisotropic dark fluid 

p = 3J(m - l)k 2 c?e mkt - l -f(R) + ^L^ +m)kt (4-27) 

which involves f(R) function. Thus if we interpret this function as a simplest 
model corresponding to ACDM, we may observe that at the initial epoch, 
the energy density is proportional to a constant vacuum energy density A 
only. However, with the passing time, it acquires dynamical terms describing 
matter densities of fluids filling the universe. These terms will expand to 
infinity in the future evolution of the universe. This indicates that the energy 
density is increasing in the future due to creation of some phantom DE and 
decay of some other components of energy in the universe, like, cold dark 
matter. Consequently, expansion of the universe is accelerated forever. 

Inserting the values of A, F, p in Eq. fl4.13p . we obtain skewness parameter 

S 

5 = tfELLEl ^— . (4.28) 

4[18Pk 2 {m - l)4 +m e mkt - 4ciq (2+m) e-( 6+m ) fc * - 3lcff{R)} 

Using this value of 5 alongwith scale factors and p in Eqs . (14 . 1 7|) - (14 . 1 9 j) and 
solving them simultaneously, the anisotropic EoS parameter u will become 

6(1 - m)c\ +m l 2 k 2 e mkt - 8c 5 lmke~ 3kt + 6lc 3 J(R) 
OJ — 77T- — > . (4.29) 

18l 2 k 2 {m - l)cf +m e mkt - 4c|cr ( +m) e-( 6+m ) fct - 3lcff{R) 

Since both the parameters contain p and hence f(R) function, it may not be 
possible to discuss their behavior independently in the future evolution. The 
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scalar curvature for exponential solution becomes 



R = -10fc 2 - - 




1 



e 



2(3+m)kt 



+ 




2a 2 -2fc*+f 



4 e 5 
3 I(3+m)fcc' 




(4.30) 



This procedure can be extended to the solution of the field equations for 
power-law volumetric expansion. 

5 Outlook 

This paper is devoted to study the exact solutions of the Bianchi type VI 
universe in the metric f(R) gravity and to discuss the recent cosmic ac- 
celeration. Two types of non-vacuum solutions are found corresponding to 
isotropic and anisotropic fluids respectively. The physical behavior of these 
solutions is discussed at early and late times of the universe, which can be 
summarized as follows: 

• The scale factors in the isotropic fluid case are zero at t = which 
shows that the spacetime exhibits point type singularity and continues 
to expand till t — > oo. In anisotropic case, the scale factors admit 
constant value at initial epoch and then start increasing with cosmic 
time approaching to very large values as t — > oo. Thus, the spacetime 
does not show any type of initial singularity in this case. 

• The model starts with physical parameters such that all being infinite 
at early times and approaches to zero at later times in the perfect 
fluid case. In the anisotropic case, the expansion scalar and the mean 
Hubble parameter are constant indicating homogenous expansion of 
the universe. The volume of the universe will increase with time due 
to expansion for both the solutions. 

• The matter density and pressure of the fluid are related by the EoS, 
where EoS parameter is used to characterize the DE into different ex- 
pansion histories. It is found that in isotropic solution, the DE has 
large negative pressure with u < — 1 which shows that the universe 
passes through phantom region. It is mentioned here that our con- 
clusions support the observational evidence of a recent supernova data 
[53]- [ES]- However, the energy density of the anisotropic fluids involves 
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f{R) function whose interpretation is given in the context of AC DM. 
In this case, the EoS parameter turns out to be time dependent. 



• In the anisotropic solution, the deviation from the isotropy along z- 
axis is zero which shows that pressure of DE along y and z axes are 
same. The skewness parameter along x-axis does not vanish even at 
the later times. However, for m > —3, the anisotropy parameter of 
the expansion becomes zero at infinite times showing that the universe 
expands isotropically regardless of the anisotropy in the fluid. It is 
mentioned here that in GR, the anisotropic fluid as well as the model 
attain isotropy in the future evolution of the universe [58j [7T] but in 
f{R) theory, the model and the fluid remain anisotropic. Further, in 
the isotropic fluid case, the universe does not achieve isotropy and the 
model expands continuously with non-zero shear scalar. It is remarked 
here that the above mentioned aspect of the model is compatible with 
existing results of GR [60, 62J. 
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